package top.jfunc.common.utils;

import org.gavaghan.geodesy.Ellipsoid;
import org.gavaghan.geodesy.GeodeticCalculator;
import org.gavaghan.geodesy.GeodeticCurve;
import org.gavaghan.geodesy.GlobalCoordinates;

import java.awt.geom.Point2D;
import java.util.List;

/**
 * 判断一个点是否在一个多边形内
 * 提供了三种计算方式，孰好孰坏在实践中比较一下
 * @author xiongshiyan at 2019/1/11 , contact me with email yanshixiong@126.com or phone 15208384257
 */
public class GeometryUtil {
    private static final int POLYGON_MIN_SIZE = 3;

    public static boolean isPtInPolygon (Point2D.Double point , List<Point2D.Double> polygon) {
        assertParams(point, polygon);

        int iSum,iIndex;
        double dLon1 , dLon2 , dLat1 , dLat2 , dLon;
        int size = polygon.size();
        iSum = 0;
        for (iIndex = 0; iIndex<size; iIndex++) {
            if (iIndex == size - 1) {
                dLon1 = polygon.get(iIndex).getX();
                dLat1 = polygon.get(iIndex).getY();
                dLon2 = polygon.get(0).getX();
                dLat2 = polygon.get(0).getY();
            } else {
                dLon1 = polygon.get(iIndex).getX();
                dLat1 = polygon.get(iIndex).getY();
                dLon2 = polygon.get(iIndex + 1).getX();
                dLat2 = polygon.get(iIndex + 1).getY();
            }
            // 以下语句判断A点是否在边的两端点的水平平行线之间，在则可能有交点，开始判断交点是否在左射线上
            if (((point.y >= dLat1) && (point.y < dLat2))
                    || ((point.y >= dLat2) && (point.y < dLat1))) {
                if (Math.abs(dLat1 - dLat2) > 0) {
                    //得到 A点向左射线与边的交点的x坐标：
                    dLon = dLon1 - ((dLon1 - dLon2) * (dLat1 - point.y) ) / (dLat1 - dLat2);
                    // 如果交点在A点左侧（说明是做射线与 边的交点），则射线与边的全部交点数加一：
                    if (dLon < point.x) {
                        iSum++;
                    }
                }
            }
        }
        return (iSum % 2) != 0;
    }

    /**
     * 返回一个点是否在一个多边形区域内， 如果点位于多边形的顶点或边上，不算做点在多边形内，返回false
     */
    public static boolean isPointInPoly(Point2D.Double point, List<Point2D.Double> polygon) {
        assertParams(point, polygon);

        java.awt.geom.GeneralPath p = new java.awt.geom.GeneralPath();
        Point2D.Double first = polygon.get(0);
        p.moveTo(first.x, first.y);
        int size = polygon.size();
        for (int i = 1; i < size; i++) {
            Point2D.Double pa = polygon.get(i);
            p.lineTo(pa.x, pa.y);
        }
        p.lineTo(first.x, first.y);
        p.closePath();
        return p.contains(point);
    }

    /**
     * 判断点是否在多边形内，如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回true
     * @param point 检测点
     * @param polygon   多边形的顶点
     * @return      点在多边形内返回true,否则返回false
     */
    public static boolean isPtInPoly(Point2D.Double point, List<Point2D.Double> polygon){
        assertParams(point, polygon);

        int N = polygon.size();
        //如果点位于多边形的顶点或边上，也算做点在多边形内，直接返回true

        //cross points count of x
        int intersectCount = 0;
        //浮点类型计算时候与0比较时候的容差
        double precision = 2e-10;
        //neighbour bound vertices
        Point2D.Double p1, p2;
        //当前点
        Point2D.Double p = point;

        //left vertex
        p1 = polygon.get(0);
        //check all rays
        for(int i = 1; i <= N; ++i){
            if(p.equals(p1)){
                //p is an vertex
                return true;
            }

            //right vertex
            p2 = polygon.get(i % N);
            //ray is outside of our interests
            if(p.x < Math.min(p1.x, p2.x) || p.x > Math.max(p1.x, p2.x)){
                p1 = p2;
                //next ray left point
                continue;
            }

            //ray is crossing over by the algorithm (common part of)
            if(p.x > Math.min(p1.x, p2.x) && p.x < Math.max(p1.x, p2.x)){
                //x is before of ray
                if(p.y <= Math.max(p1.y, p2.y)){
                    //overlies on a horizontal ray
                    if(p1.x == p2.x && p.y >= Math.min(p1.y, p2.y)){
                        return true;
                    }

                    //ray is vertical
                    if(p1.y == p2.y){
                        //overlies on a vertical ray
                        if(p1.y == p.y){
                            return true;
                            //before ray
                        }else{
                            ++intersectCount;
                        }
                        //cross point on the left side
                    }else{
                        //cross point of y
                        double xinters = (p.x - p1.x) * (p2.y - p1.y) / (p2.x - p1.x) + p1.y;
                        //overlies on a ray
                        if(Math.abs(p.y - xinters) < precision){
                            return true;
                        }

                        //before ray
                        if(p.y < xinters){
                            ++intersectCount;
                        }
                    }
                }
                //special case when ray is crossing through the vertex
            }else{
                //p crossing over p2
                if(p.x == p2.x && p.y <= p2.y){
                    //next vertex
                    Point2D.Double p3 = polygon.get((i+1) % N);
                    //p.x lies between p1.x & p3.x
                    if(p.x >= Math.min(p1.x, p3.x) && p.x <= Math.max(p1.x, p3.x)){
                        ++intersectCount;
                    }else{
                        intersectCount += 2;
                    }
                }
            }
            //next ray left point
            p1 = p2;
        }

        //偶数在多边形外
        return intersectCount % 2 != 0;
    }
    private static void assertParams(Point2D.Double point, List<Point2D.Double> polygon) {
        if(null == point || null == polygon || polygon.size() < POLYGON_MIN_SIZE) {
            throw new IllegalArgumentException("参数不能为空，且多边形点数大于3");
        }
    }



    /**
     * 默认地球半径
     */
    private static final double EARTH_RADIUS = 6_378_137D; // 赤道半径(单位m)

    /**
     * 功能描述：计算两个经纬度的距离
     * 基于googleMap中的算法得到两经纬度之间的距离, 计算精度与谷歌地图的距离精度差不多
     * @param longitude1 第一点的精度
     * @param latitude1 第一点的纬度
     * @param longitude2 第二点的精度
     * @param latitude2 第二点的纬度
     * @return 返回的距离，单位m
     * @author cakin
     * @serialData  2021/10/30
     */
    public static double googleDistance(double longitude1, double latitude1, double longitude2, double latitude2) {
        double radLat1 = Math.toRadians(latitude1);
        double radLat2 = Math.toRadians(latitude2);
        double a = radLat1 - radLat2;
        double b = Math.toRadians(longitude1) - Math.toRadians(longitude2);
        double s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2) + Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin(b / 2), 2)));
        return Math.round(s * EARTH_RADIUS * 10_000D) / 10_000D;
    }

    /**
     * （反余弦计算方式）
     *
     * @param longitude1 第一点的经度
     * @param latitude1  第一点的纬度
     * @param longitude2 第二点的经度
     * @param latitude2  第二点的纬度
     * @return 返回的距离，单位m
     */
    public static double acosDistance(double longitude1, double latitude1, double longitude2, double latitude2) {
        // 经纬度（角度）转弧度。弧度作为作参数，用以调用Math.cos和Math.sin
        // A经弧度
        double radiansAX = Math.toRadians(longitude1);
        // A纬弧度
        double radiansAY = Math.toRadians(latitude1);
        // B经弧度
        double radiansBX = Math.toRadians(longitude2);
        // B纬弧度
        double radiansBY = Math.toRadians(latitude2);

        // 公式中“cosβ1cosβ2cos（α1-α2）+sinβ1sinβ2”的部分，得到∠AOB的cos值
        double cos = Math.cos(radiansAY) * Math.cos(radiansBY) * Math.cos(radiansAX - radiansBX) + Math.sin(radiansAY) * Math.sin(radiansBY);

        // 反余弦值
        double acos = Math.acos(cos);

        // 最终结果
        return EARTH_RADIUS * acos;
    }


    /**
     * （高德地图计算方法）
     *
     * @param longitude1 第一点的经度
     * @param latitude1  第一点的纬度
     * @param longitude2 第二点的经度
     * @param latitude2  第二点的纬度
     * @return 返回的距离，单位m
     */
    public static Double amapDistance(double longitude1, double latitude1, double longitude2, double latitude2) {
        longitude1 *= 0.01745329251994329;
        latitude1 *= 0.01745329251994329;
        longitude2 *= 0.01745329251994329;
        latitude2 *= 0.01745329251994329;
        double var1 = Math.sin(longitude1);
        double var2 = Math.sin(latitude1);
        double var3 = Math.cos(longitude1);
        double var4 = Math.cos(latitude1);
        double var5 = Math.sin(longitude2);
        double var6 = Math.sin(latitude2);
        double var7 = Math.cos(longitude2);
        double var8 = Math.cos(latitude2);
        double[] var10 = new double[3];
        double[] var20 = new double[3];
        var10[0] = var4 * var3;
        var10[1] = var4 * var1;
        var10[2] = var2;
        var20[0] = var8 * var7;
        var20[1] = var8 * var5;
        var20[2] = var6;

        return Math.asin(Math.sqrt((var10[0] - var20[0]) * (var10[0] - var20[0]) + (var10[1] - var20[1]) * (var10[1] - var20[1]) + (var10[2] - var20[2]) * (var10[2] - var20[2])) / 2.0) * 1.27420015798544E7;
    }


    /**
     * （利用第三方jar包计算）
     * 计算两个经纬度之间的距离
     *https://mvnrepository.com/artifact/org.gavaghan/geodesy
     *compile 'org.gavaghan:geodesy:1.1.3'
     * @param longitude1 第一点的经度
     * @param latitude1  第一点的纬度
     * @param longitude2 第二点的经度
     * @param latitude2  第二点的纬度
     * @param ellipsoid  计算方式
     * @return 返回的距离，单位m
     */
    public static double getDistance(double longitude1, double latitude1, double longitude2, double latitude2, Ellipsoid ellipsoid) {
        // 创建GeodeticCalculator，调用计算方法，传入坐标系、经纬度用于计算距离
        GeodeticCurve geoCurve = new GeodeticCalculator().calculateGeodeticCurve(ellipsoid,
                new GlobalCoordinates(latitude1, longitude1),
                new GlobalCoordinates(latitude2, longitude2));
        return geoCurve.getEllipsoidalDistance();
    }

    public final static double x_pi = Math.PI * 3000.0 / 180.0;
    /**
     * 高德坐标转百度坐标
     * - 高德地图使用的是GCJ-02坐标系（火星坐标系）
     * - 百度地图使用的是BD-09坐标系
     */
    public static LngLatPoint gd2bd(LngLatPoint gdPoint) {
        double x = gdPoint.getLongitude(), y = gdPoint.getLatitude();
        double z = Math.sqrt(x * x + y * y) + 0.00002 * Math.sin(y * x_pi);
        double theta = Math.atan2(y, x) + 0.000003 * Math.cos(x * x_pi);
        double bd_lon = z * Math.cos(theta) + 0.0065;
        double bd_lat = z * Math.sin(theta) + 0.006;
        return new LngLatPoint(bd_lon, bd_lat);
    }

    /**
     * 百度坐标转换高德坐标
     * - 高德地图使用的是GCJ-02坐标系（火星坐标系）
     * - 百度地图使用的是BD-09坐标系
     */
    public static LngLatPoint bd2gd(LngLatPoint bdPoint) {
        double x = bdPoint.getLongitude() - 0.0065, y = bdPoint.getLatitude() - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * x_pi);
        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * x_pi);
        double gd_lon = z * Math.cos(theta);
        double gd_lat = z * Math.sin(theta);
        return new LngLatPoint(gd_lon, gd_lat);
    }
}